The present invention relates to a method to migrate seismic time/amplitude data into depth/amplitude data. As explained by Robinson and Treitel in geophysical signal analysis, at page 374, the process by which the apparent dips in a seismic record section are converted to time dips is known as migration. As explained on page 378 of geophysical signal analysis, migration is the process of constructing the reflector surface (the actual interface of an underground geologic structure that produces a reflection) from the record surface (the surface defined by reflections before accounting for slope). More specifically, depth migration may also be viewed as the process of placing reflection events recorded in the seismic data at their proper spatial location. This is a computationally intensive process especially when the velocity field varies laterally at any given depth. The present invention relates to improving the computation speed of depth migration using wavelet analysis.
The term "downward continuation" is a term known in the art. Downward continuation is the computation of the wavefield, knowing its value at a known depth (e.g., depth=z=0), down to any given depth, using the wave equation. That is, downward continuation is the computation of the wavefield from where it is known, on the surface, to where it has not been recorded, at some depth within the earth, assuming some velocity model. In the following description of the present invention, we have chosen to use the explicit formulation of the downward continuation operator from depth z to depth z+.DELTA.z, expressed in the wavelet space.
The seismic imaging art includes a variety of techniques for improving definition in the seismic image. For example, Whitmore, Jr. et al., U.S. Pat. No. 4,766,574, relates to a method for imaging multi-component seismic data to obtain better depth images of the earth's subsurface geological structure as well as better estimates of compressional and shear wave interval velocities. FIG. 10 of Whitmore, Jr. et al. depicts a generalized process flow diagram. The present invention provides a new method for processing seismic data that would find application in step 50 in FIG. 10. In Whitmore, Jr. et al., the combination of equation (3), the imaging principle, and equations (4) and (5), wavefield continuation, constitutes depth migration, as undertaken at step 50 of FIG. 10. The incident and reflected wavefield are mathematically continued down (i.e., the downward continuation problem) into the earth's subsurface via equations (4) and (5) and a reflectivity is computed from equation (3). The migration process as contemplated by Whitmore, Jr. et al. and other known migration techniques is extremely calculation intensive when the velocity field varies laterally at a given depth and is consequently time consuming. It should be noted that the present invention is generally applicable to depth migration, independent of Whitmore, Jr. et al., but that Whitmore, Jr. et al. would benefit significantly from the application of the present invention to their step 50.
The present invention utilizes wavelet transforms and the wavelet space to perform the calculations of the downward continuation (or migration) process. Wavelet analysis is a relatively recently developed field of mathematical and signal analysis. The literature is now abundant in applications of wavelet transforms in the field of image processing, speech processing, and telecommunications. In contrast, the literature is sparse in geophysics.
Geophysicists generally use the term "wavelet" to refer to the waveform containing the seismic source characteristics. In the context of wavelet analysis, the word "wavelet" refers to families of mathematical functions obtained by translation (usually by an integer) and dilation (usually by powers of two) of a single function. (As used herein, the term "translation" refers to the movement of the center time of the wavelet function and the terms "dilation" and "contraction" refer to the adjustment of the width of the wavelet.) The following discussion introduces the scaling function and wavelets of wavelet analysis.
A family of functions .PHI..sub.nk, whose integral is unity, is generated from a single function .PHI., by dilation or contraction by powers of two as well as translation, through the relation: EQU .PHI..sub.nk (x)=2.sup.-n/2 .PHI.(x/2.sup.n -k) (1)
.PHI. is called a scaling function. The scaling function satisfies the relation: EQU .PHI.(x)=.SIGMA..sub.k .alpha..sub.k .PHI.(2x-k) (2)
.PHI. is compactly supported (as that term is known from Daubechies, Orthonormal Bases of Compactly Supported Wavelets, Communication of Pure and Applied Mathematics, Vol. XLI, pp. 909-996), or identically zero outside of a finite length interval when a finite number of a.sub.k coefficients is non-zero. The wavelet is an associated function whose integral is zero. It is obtained through: EQU .PSI.(x)=.SIGMA..sub.k (-1).sup.k .alpha..sub.1-k .PHI.(2x-k) (3)
In a manner similar to the family of functions defined in equation (1), .PSI..sub.nk functions are generated from .PSI. by dilation and translation. It has been demonstrated that, for some choices of the a.sub.k coefficients, the .PSI..sub.nk functions form an orthonormal basis for square integrahie functions on the real axis. A discrete wavelet transform is the mapping of any square integrahie function f into its wavelet coefficients: EQU f(x)=.SIGMA..sub.nk (f,.PSI..sub.nk).PSI..sub.nk (x) (4) where: EQU f.sub.nk =(f,.PSI..sub.nk)=.intg.f(x).PSI..sub.nk (x)dx (5)
are the wavelet coefficients. The integers n and k span the whole range of integers.
Thus, the word "wavelet" refers to families of mathematical functions obtained by translation and dilation of a single function, possessing in addition the property of generating a basis for the space of square integrable functions. For further explanation of wavelet transforms, consider Lecture Notes from the American Mathematical Society Short Course Series on Wavelet and Applications, Jan. 11-12, 1993, San Antonio, Tex.
As described in Wavelets and Signal Processing, O. Rioul and M. Vetterli, IEEE SP Magazine, October, 1991, the wavelet transform (WT) provides an alternative to the short-time Fourier transform (STFT). In contrast to the STFT, which uses a single analysis window, the WT uses short windows at high frequencies and long windows at low frequencies (i.e., dilations and contractions), thus providing constant relative bandwidth frequency analysis. It may be desirable to see the WT as a signal decomposition onto a set of basis functions. They are obtained from a single prototype wavelet by dilations and contractions (scaling) as well as shifts or translations.
The set of coefficients corresponding to one given level n (Equation 5) can be interpreted as a projection of f onto a subspace W.sub.n of square integrable functions. (As used herein, the term "level" will best be understood with reference to FIG. 5 and the discussion below with regard to that figure.) W.sub.n is generated by the .PSI..sub.nk functions. The related subspace V.sub.n is generated by the associated .PSI..sub.nk functions. The two subspaces are disjoint and their direct sum (or union) is V.sub.n-1, the subspace at the finer scale. The projection of f onto V.sub.n where n is decreasing may be intuitively interpreted as a finer and finer approximation of f.
Thus, the projection of a function at a given level may be decomposed into a smoothly varying part containing its broad features (i.e., the low frequency part), and a high frequency part with its sharp features (or edges).
Daubechies, in Orthonormal Bases of Compactly Supported Wavelets, Communication of Pure and Applied Mathematics, Vol. XLI, pp. 909-996, studies a family of such wavelets with the properties of compact support (.PHI. and .PSI. vanish outside of a finite interval) and orthonormality. The computation of the coefficients h.sub.nk, which constitute the wavelet transform, may be obtained efficiently, without evaluating the integral in equation (5), above. The cost of the transform (i.e., the computation overhead) is proportional to N, where N is the number of points of the input data.
Thus, the depth migration process being computationally intensive, there is a need to reduce the number of data points that must be analyzed to provide an adequate seismic image to thereby reduce time and cost involved. Wavelet analysis provides this reduction in the number of data points needed in the analysis.